%% SCRIPT TO MERGE MARSHALL FIRE DOW RADAR AND WRF METEOROLOGICAL CROSS SECTION DATA

%produced by Neil P. Lareau summer 2022
%for questions contact: nlareau@unr.edu

close all
clear all

%% LOAD SOME OF THE RADAR DATA
load('marhsall_plume_d1.mat'); %this can be 'd1' or 'd2' for deployment 1 or 2
xv=unique(xxx(:)); %extract unique x-distance values
zv=unique(zzz(:)); %extract unique y-distance values
alt=1 %place holder
lons=stationlon; %station longitude 
lats=stationlat; %station latitude
re=6.37*10^6; %earth radius
scalefactor=(re.*cosd(nanmean(lats))).*(pi./180); %scale factor for converting x and longitude data

%% LOAD WRF DATA 

fname='/Users/nlareau/Dropbox/MARSHALL_FIRE/wrfout_d02_2021-12-30_22:15:00'; %file path to wrf output file
ncdisp(fname) %display netcdf contents
wtime=ncread(fname,'Times'); %times
T=ncread(fname,'T'); %temperature


%% examine surface map
latsW=squeeze(double(ncread(fname,'XLAT'))); %wrf latitude data
lonsW=squeeze(double(ncread(fname,'XLONG'))); %wrf longitude data

thgt=squeeze(double(ncread(fname,'HGT'))); %terrain heights
%visualize WRF terrain
figure(1);clf;
pcolor(lonsW,latsW,thgt);shading flat;
hold on;
plot(stationlon,stationlat,'*k'); %add radar location
plot(lonsW(:,300),latsW(:,300),'--k') %plot an east-west line in the region of interest

%% Contruct geopotential data
gph=(squeeze(double(ncread(fname,'PH')))+squeeze(double(ncread(fname,'PHB'))))./9.81; %base state plus perturbation geopotential 

% produce a quick cross section
figure(2);clf;
xcross=repmat(lonsW(:,300),1,size(T,3)); %longitude slice
tcross=squeeze(nanmean(T(:,280:320,:),2)); %temperatreu slice (mean over a subset of meridional points spanning the radar location)
hgtcross=squeeze(nanmean(gph(:,280:320,1:end-1),2)); %geopotential heights over the same region
pcolor(xcross,hgtcross,tcross+300);shading flat; %temperature cross section in Kelvin
ylim([1200 7000])
caxis([290 330]);
hold on;
contour(xcross,hgtcross,tcross+300,[270:2:500],'k')
cbh=colorbar;

%% add wind cross section
uwind=double(squeeze(ncread(fname,'U')));
addpath('~/Dropbox/MATLAB/SGP_LIDAR_CODES/')
figure(3);clf;
xcross=repmat(lonsW(:,300),1,size(T,3));
ucross=squeeze(nanmean(uwind(:,300,:),2));
hgtcross=squeeze(nanmean(gph(:,300,1:end-1),2));
pcolor(xcross,hgtcross,ucross(1:end-1,:));shading flat;
ylim([1200 7000])
caxis([-2 40]);
hold on;
contour(xcross,hgtcross,tcross+300,[270:2:500],'k')
colormap(rbmapper(40,-2))


%% Detailed cros section for publication
[~,tidx]=min(abs(timemaster-datenum(2021,12,30,22,17,30)));

for ttt=tidx %was set up as a loop but now just uses the target time defined above in tidx, can be revereted if you want

    addpath('~/Dropbox/MATLAB/')
    load('radar_cmap.mat') %radar colormap
    figure(5);clf; set(gcf,'color','w','pos',[100 100 1000 1000]) %call figure
    subplot(3,1,1:2) %upper subplot planels
    dbzbar=squeeze(nanmax(nanmax(dbzmaster((ttt-4):(ttt+4),:,:,:),[],2),[],1));%meridional and time maximum of the radar reflectivity
    dbzbar(dbzbar<-10)=nan; %remove low values
    dbzx=nanmax(dbzbar(:,5:end),[],2) %this is the "column" maximum from the time-meridional maximum above
    pcolor((xv./scalefactor)+stationlon,zv+nanmean(stationalt),medfilt2(dbzbar)'); %display the time-meridional maximum and smooth with 3x3 median filter
    shading flat;
    colormap(cmap)
    caxis([-20 35]);
    cbh=colorbar; ylabel(cbh,'dbZe','fontsize',15,'fontweight','bold')
    grid on; box on; set(gca,'layer','top','linewidth',2,'fontsize',18,'fontweight','bold')
    xlabel('Longitude')
    ylabel('Alt [m MSL]')
    pause(.1)
    xlim([-105.42 -104.8])
    hold on;
    ylim([1400 6000])
    text(-105.4,5600,'a)','fontsize',20,'fontweight','bold')
    title(strcat(datestr(timemaster(ttt-4),'HHMM'),'-',datestr(timemaster(ttt+4),'HHMM'),{' UTC'}));

    %add the terrain cross section ontop of the figure
    ah=area(xcross(:,1),thgt(:,300));
    ah.LineWidth=3;%,'k','linewidth',3)
    ah.FaceColor=[.7 .7 .7];
    hold on;

    %% add potential temperature contuors at 1 and 5 K intervals. Label the 5 K values
    contour(xcross,hgtcross,tcross+300,[270:1:500],'k')
    [c,h]=contour(xcross,hgtcross,tcross+300,[270:5:500],'k','linewidth',2);
    clabel(c,h,'LabelSpacing',2000,'fontsize',20,'fontweight','bold')

    %lower subplot showing along wind change in column maximum reflectivity
    subplot(3,1,3);cla;hold on;
    dbzxc=dbzx;
    dbzxc(isnan(dbzxc))=-30;%fill nan values before smooth
    [bb,cc]=butter(5,1/5); %by direction 5th order butterworth filter
    dbzxsm=filtfilt(bb,cc,dbzxc); %smooth it... this shouldn't shift peaks
    dbzxsm(dbzxsm<0)=nan; %replace nans
    plot((xv./scalefactor)+stationlon,dbzx./1) %plot the "raw"
    plot((xv./scalefactor)+stationlon,dbzxsm./1,'linewidth',4) %plot the smoothed
    xlim([-105.42 -104.8]) %xlimits
    cbh=colorbar; cbh.Visible='off';
    grid on; box on; set(gca,'layer','top','linewidth',2,'fontsize',18,'fontweight','bold')
    xlabel('Longitude')
    ylabel('dbZe')
    ylim([-2 32]./1)
    text(-105.4,29,'b)','fontsize',20,'fontweight','bold')

end


%% Now produce a similar spectrum width cross section
dbz1=dbzmaster; %store dbz in dbz1
load('marhsall_plume_d1_SW.mat') %this loads the post-processed spectrum width data... apparently I was lazy and accidentally save the spectrum width (SW) data as "dbzmaster"... hereafter dbzmaster contains the spectrum width

%% cros section
for ttt=tidx

    addpath('~/Dropbox/MATLAB/')
    load('radar_cmap.mat')
    figure(5);clf; set(gcf,'color','w','pos',[100 100 1000 1000])
    subplot(3,1,1:2)
    dbzbar=squeeze(nanmax(nanmax(dbzmaster((ttt-4):(ttt+4),:,:,:),[],2),[],1)); %this is the time-meridional maxima of the spectrum width, don't be fooled by the variable name
    dbzbar(dbzbar<1)=nan;
    dbzx=nanmax(dbzbar(:,5:end),[],2)
    pcolor((xv./scalefactor)+stationlon,zv+nanmean(stationalt),medfilt2(dbzbar)');
    shading flat;
    colormap(cmap)
    colormap(turbo(16))
    caxis([0 9]);
    cbh=colorbar; ylabel(cbh,'spectrum width [m s^{-1}]','fontsize',15,'fontweight','bold')
    grid on; box on; set(gca,'layer','top','linewidth',2,'fontsize',18,'fontweight','bold')
    xlabel('Longitude')
    ylabel('Alt [m MSL]')
    pause(.1)
    xlim([-105.42 -104.8])
    hold on;

    ylim([1400 6000])
    text(-105.4,5600,'a)','fontsize',20,'fontweight','bold')
    title(strcat(datestr(timemaster(ttt-4),'HHMM'),'-',datestr(timemaster(ttt+4),'HHMM'),{' UTC'}));
    ah=area(xcross(:,1),thgt(:,300));
    ah.LineWidth=3;%,'k','linewidth',3)
    ah.FaceColor=[.7 .7 .7];
    hold on;
    contour(xcross,hgtcross,tcross+300,[270:1:500],'k')
    [c,h]=contour(xcross,hgtcross,tcross+300,[270:5:500],'k','linewidth',2);
    clabel(c,h,'LabelSpacing',2000,'fontsize',20,'fontweight','bold')
    % colormap(rbmapper(40,-2))

    subplot(3,1,3);cla;hold on;
    dbzxc=dbzx;
    dbzxc(isnan(dbzxc))=-30;
    [bb,cc]=butter(5,1/5);
    dbzxsm=filtfilt(bb,cc,dbzxc);
    dbzxsm(dbzxsm<0)=nan;
    plot((xv./scalefactor)+stationlon,dbzx./1)
    plot((xv./scalefactor)+stationlon,dbzxsm./1,'linewidth',4)
    xlim([-105.42 -104.8])
    cbh=colorbar; cbh.Visible='off';
    grid on; box on; set(gca,'layer','top','linewidth',2,'fontsize',18,'fontweight','bold')
    xlabel('Longitude')
    ylabel('spectrum width [m s^{-1}]')
    ylim([2 10]./1)
    text(-105.4,9,'b)','fontsize',20,'fontweight','bold')

end
